%BOOTLEG PLUME VISUALIZATION AND COORDINATE ROTATION

close all
clear all


%% Constants
re=6.37*10^6;
daysec=24*60*60;


%% Load the 3D gridded radar data
% load('~/Dropbox/BOOTLEG_FIRE/bootleg_plume.mat','dbzmaster','xx','yy','zz','timer','StationLatitude','StationLongitude','StationElevationInMeters')
load('~/Dropbox/BEAR_FIRE/Bear_Fire_Plume_new.mat','dbzmaster','XX','YY','ZZ','timer','StationLatitude','StationLongitude','StationElevationInMeters')
re=6.37*10^6;
scalefactor=(re.*cosd(StationLatitude)).*(pi./180);
%% LOAD SRTM TERRAIN DATA:
addpath('~/Dropbox/MATLAB/')

X1=readhgt('~/Dropbox/BEAR_FIRE/N39W122.hgt');
[ltmt1,lnmt1]=meshgrid(X1.lat,X1.lon);
X2=readhgt('~/Dropbox/BEAR_FIRE/N39W121.hgt');
[ltmt2,lnmt2]=meshgrid(X2.lat,X2.lon);
%%
figure(1);clf; clf; hold on; pcolor(lnmt1,ltmt1,X1.z');shading flat;
pcolor(lnmt2,ltmt2,X2.z');shading flat;
%%
ltall=cat(1,ltmt1,ltmt2);
lnall=cat(1,lnmt1,lnmt2);
Zall=double(cat(1,X1.z',X2.z'));



%%
% X3=readhgt('~/Dropbox/BOOTLEG_FIRE/N42W121.hgt');
% Z3=double(X3.z');%need to rotate these data for some reason...
% [ltmt3,lnmt3]=meshgrid(X3.lat,X3.lon);  %make a lat/long mesh
xterrain=(lnall-StationLongitude).*scalefactor; %convert to x from radar site (probably KDAX)
yterrain=(ltall-StationLatitude).*111180; %convert to y from radar site (probably KDAX)
Z3=Zall;
%create terrain scattered interpolant for use in mapping x-y fire points to
%vertical points
tic
FT=TriScatteredInterp(xterrain(:),yterrain(:),Z3(:)); %FT is an interpolant into which we can pass x,y points later
toc


%% Produce 3D plume volume over terrain map
[val,tidx]=min(abs(timer-datenum(2020,9,9,2,40,0)));

for ff=tidx%68:84%54:length(timer)-1;

    f=figure(6);clf        % pcolor(terrain_lon,terrain_lat,Z);shading flat;
    set(f,'color','w','pos',[361   169   1300   770]);
    hold on;
    strd=2; %makes plotting faster, but can be set to 1 to get full terrain resolution if needed
    surf(xterrain(1:strd:end,1:strd:end),yterrain(1:strd:end,1:strd:end),Z3(1:strd:end,1:strd:end));shading flat;
    %     surf(xterrain2(1:strd:end,1:strd:end),yterrain2(1:strd:end,1:strd:end),zall2(1:strd:end,1:strd:end));shading flat;

    daspect([1 1 .6]);
    colormap([.95 1 .95]);%demcmap([0 2000]));
    caxis([00 2000]);
    hold on;
    material dull;
    grid on; box on; set(gca,'linewidth',2,'layer','top');
    camlight(40,10,'infinite')
    lighting flat;
    hold on;
    grid on; box on; set(gca,'linewidth',2,'layer','top');
    lighting gouraud;
    zlabel('Altitude [m]','fontsize',15,'fontweight','bold');
    xlabel('x-dist [m]');
    ylabel('y-dist [m]');
    set(gca,'fontsize',15,'fontweight','bold')
    %extract the current time step dbz volume
    dbzvol=squeeze(nanmax(dbzmaster(ff,:,:,:),[],1));

    dbzlist=[00 10 20 27 30 34 38]; %list of dbz surface values
    clrs=[.5 .48 .4; .6 .58 .5; 1 .7 .3; 1 .7 .3; 1 0 0; .9 0 0; .8 0 0]; %list of face colors for plumes
    trans=[.1 .2 .3 .3 .3 .5 .7]; % surface transparency

    %create transparent isosurfaces at the dbz list values
    for dbs=1:length(dbzlist)
        if dbzlist(dbs)<20
            p = patch(isosurface(XX,YY,ZZ,smooth3(dbzvol,'box',[7 7 7]),dbzlist(dbs)));
            isonormals(XX,YY,ZZ,smooth3(dbzvol,'box',[7 7 7]),p);
        else
            p = patch(isosurface(XX,YY,ZZ,smooth3(dbzvol,'box',[3 3 3]),dbzlist(dbs)));
            isonormals(XX,YY,ZZ,smooth3(dbzvol,'box',[3 3 3]),p);
        end
        p.FaceColor = clrs(dbs,:);
        p.FaceAlpha=trans(dbs);
        p.EdgeColor = 'none';
    end

    lighting gouraud; camlight(10,-100);
    zlim([0 11000]);
    daspect([1 1 .5]);
    xlim([0 8]*10^4);
    ylim([.7 1.6]*10^5);
    title(datestr(nanmean(timer(ff)),'HH:MM'));
    view(120,12);
end

%% Plan view map
dbzmap=squeeze(nanmax(nanmax(dbzmaster(ff-2:ff+2,:,:,5:end),[],1),[],4));
figure(74);clf;set(gcf,'color','w','pos',[100 100 1000 1000])
dbzmap(dbzmap<0)=nan;
pcolor(squeeze(XX(:,:,1)), squeeze(YY(:,:,1)),medfilt2(dbzmap,[5 5]));shading flat; caxis([0 40])
hold on;
daspect([1 1 1])

% add line showing the spot distnaces
[x,y]=ginput(2)
plot(x,y,'--k','linewidth',4)
dister=sqrt((x(2)-x(1)).^2 + (y(2)-y(1)).^2);
theta=atand((y(2)-y(1))./(x(2)-x(1)));
xticks=get(gca,'XTick');set(gca,'XTickLabel',xticks/1000)
yticks=get(gca,'YTick');set(gca,'YTickLabel',yticks/1000)
xlabel('distance [km]')
ylabel('distance [km]')
grid on; box on;
set(gca,'layer','top','linewidth',2,'fontsize',18,'fontweight','bold')
colormap(turbo(15));%[flipud(bone(48)); hot(15)])
cbh=colorbar;ylabel(cbh,'dbZ');
caxis([10 40]);
hold on;
% contour(xterrain,yterrain,medfilt2(Z3,[10 10]),[0:100:3000],'k')


%% axis rotation to create x distances along wind including spot fire distances

%start with initial x,y,z matrices and dbzvolume data
xvec=unique(XX(:));%create and xvector from the x matrix
yvec=unique(YY(:)); %construct the y vector from the ymatrix
zvec=unique(ZZ(:));%create the vertical vector

figure(100);clf;hold on;set(gcf,'pos',[100 100 1000 500])
subplot(1,3,1);cla; hold on;
dbzmap(isnan(dbzmap))=-30;
daspect([1 1 1])
hold on; contour(xvec,yvec,medfilt2(dbzmap),[5:5:50],'k');
% add an arbitrary point to test the rotation for spot fire data
[xsp,ysp]=ginput(1)
plot(xsp,ysp,'sr','markersize',20)
%now recenter
[x0,y0]=ginput(1)
xv2=xvec-x0;
yv2=yvec-y0;
subplot(1,3,2);cla; hold on;
daspect([1 1 1])
hold on; contour(xv2,yv2,medfilt2(dbzmap),[5:5:50],'k');


%% now rotate axes through angle theta
% theta=28;
[xvm,yvm,zvm]=meshgrid(xv2,yv2,zvec);
x2=xvm.*cosd(theta) + yvm.*sind(theta);
y2=-xvm.*sind(theta) + yvm.*cosd(theta);
subplot(1,3,3);cla; hold on;
daspect([1 1 1])
hold on; contour(squeeze(x2(:,:,1)),squeeze(y2(:,:,1)),medfilt2(dbzmap),[5:5:50],'k');

%rotate the spot fire point
xsp2=(xsp-x0).*cosd(theta) +(ysp-y0).*sind(theta);
ysp2=-(xsp-x0).*sind(theta) + (ysp-y0).*cosd(theta);
plot(xsp2,ysp2,'sr','markersize',20)
x02=x0;
y02=y0;


%% Next interpolate the rotated data to a regularly spaced x',y',z' catersian grid so we can average over the across plume direction easily

xvnew=[-4*10^4:250:4*10^4];
yvnew=[-3*10^4:250:3*10^4];
zvnew=[500:250:13000];
[xxx,yyy,zzz]=meshgrid(xvnew,yvnew,zvnew);

% tidx=find(timer>datenum(2020,9,8,19,0,0) & timer<datenum(2020,9,9,4,00,0));
tidx=find(timer>datenum(2020,9,9,2,30,0) & timer<datenum(2020,9,9,2,55,0));

dbzrot=nan(length(tidx),size(xxx,1),size(xxx,2),size(xxx,3));
timerot=nan(length(tidx),1);
for ttt=1:length(tidx)
    disp(ttt./length(tidx))
    dbzvol=squeeze(dbzmaster(tidx(ttt),:,:,:));
    tic
    %     F1=TriScatteredInterp(x2(:),y2(:),zvm(:),dbzvol(:));
    F1=TriScatteredInterp(x2(dbzvol>0),y2(dbzvol>0),zvm(dbzvol>0),dbzvol(dbzvol>0));
    dbzvint=F1(xxx(:),yyy(:),zzz(:));
    dbzvint=reshape(dbzvint,size(xxx));
    toc
    dbzrot(ttt,:,:,:)=dbzvint;
    timerot(ttt)=timer(tidx(ttt));

end
%%
dbzmap2=squeeze(nanmax(dbzvint,[],3));
hold on; contour(squeeze(xxx(:,:,1)),squeeze(yyy(:,:,1)),medfilt2(dbzmap2,[5 5]),[5:5:50],'m');

%% now rotate the terrain and remap to the same x,y grid as the radar

xtr=(xterrain-x0).*cosd(theta) + (yterrain-y0).*sind(theta);
ytr= -(xterrain-x0).*sind(theta) + (yterrain-y0).*cosd(theta);

figure(28);clf;
subplot(1,3,1);
pcolor(xterrain,yterrain,Z3); shading flat;
hold on; contour(xvec,yvec,medfilt2(dbzmap),[5:5:50],'k');

daspect([1 1 1])
subplot(1,3,2);cla;
pcolor(xtr,ytr,Z3);shading flat;
daspect([1 1 1])
hold on;
hold on; contour(squeeze(x2(:,:,1)),squeeze(y2(:,:,1)),medfilt2(dbzmap),[5:5:50],'k');

tic
F2=TriScatteredInterp(xtr(1:1:end)',ytr(1:1:end)',double(Z3(1:1:end))');
toc

xsfc=squeeze(xxx(:,:,1));
ysfc=squeeze(yyy(:,:,1));
terrainint=F2(xsfc(:),ysfc(:))
terrainint=reshape(terrainint,size(xsfc));
subplot(1,3,3);cla; hold on;
pcolor(xsfc,ysfc,terrainint);shading flat;
hold on; contour(squeeze(x2(:,:,1)),squeeze(y2(:,:,1)),medfilt2(dbzmap,[9 9]),[5:5:50],'k');
daspect([1 1 1])

%% Rotate Fire perimeter data

load('~/Dropbox/RADAR_PERIM_TRACKING/bear_fire_auto_perimeters_20211029.mat','fire_cloud');
ftimer=[fire_cloud.time];
% [val,ttidx]=min(abs(ftimer-datenum(2020,9,9,1,10,0)));
%
% %note these data are in x,y coordinates relative to KBBX, whereas the other
% %radar data is relaive to KDAX. Thus we need to add an xoffset and yoffset
% %to these data
% StationLatitudekbbx              = 39.4961;
% StationLongitudekbbx             = -121.6317;
% StationElevationInMeterskbbx    = 53;
% xoff=-(StationLongitude-StationLongitudekbbx).*scalefactor;
% yoff=-(StationLatitude-StationLatitudekbbx).*111180;
%
% xnowl=(fire_cloud(ttidx).xperim./scalefactor)+StationLongitudekbbx;
% ynowl=(fire_cloud(ttidx).yperim./111180)+StationLatitudekbbx;
%
% xnow=(xnowl-StationLongitude).*scalefactor;
% ynow=(ynowl-StationLatitude).*111180;
% xrot=((xnow-x0).*cosd(theta)) + ((ynow-y0).*sind(theta));
% yrot=-(xnow-x0).*sind(theta) + ((ynow-y0).*cosd(theta));
% fz=F2(xrot,yrot);
%
%
% % sp=scatter(spot_xr(alltimes<tideal(vvv-1)),spotz(alltimes<tideal(vvv-1)),300,'^k','filled','markeredgecolor','k')
% % sp=scatter(spot_xr(alltimes<tideal(vvv+1) & alltimes>tideal(vvv)),spotz(alltimes<tideal(vvv+1) & alltimes>tideal(vvv)),300,'^m','filled','markeredgecolor','k')
% % figure(79);clf;hold on;
% tcross=squeeze(terrainint(yvnew==0,:));
% tinner=inpolygon(xvnew,zeros(size(xvnew)),xrot,yrot);
%     sp3=scatter(xrot,fz,150,'sr','filled')
% scatter(xvnew(tinner),tcross(tinner),150,'sr','filled')
grid on; box on;
set(gca,'layer','top','linewidth',2,'fontsize',18,'fontweight','bold')

xticks=get(gca,'Xtick');yticks=get(gca,'Ytick');
set(gca,'XTickLabel',xticks/1000)
xlabel('Dist [km]');
% title(datestr(tideal(vvv)-8/24,'yyyy/mm/dd HH:MM PM'))
pause(3)


%% new cross section
close all
aa=0
axx=tight_subplot(1,1)
for ttt=3%1:1:(size(dbzrot,1))
    close all;
    figure(500);clf; set(gcf,'color','w','pos',[100 100 1500 1500])
%     subplot(1,6,1:4);cla; hold on;
    aa=aa+1;
    dbznow=squeeze(nanmax(dbzrot(ttt,:,:,:),[],1));
    dbznow(dbznow<10)=nan;
    dbznow(isnan(dbznow))=-30;
    dbzcross=medfilt2(squeeze(nanmax(dbznow,[],1))',[1 1]);
    dbzcross(dbzcross<0)=nan;
    load('radar_cmap.mat')
    pcolor(xvnew-(-2.8*10^4),zvnew,dbzcross); shading flat;
    caxis([10 48])
    load('~/Dropbox/BEECHIE_CREEK/BEACHI_SMOKE_CMAP.mat')
    colormap(BEACHIE_SMOKE);
    xlim([0 6]*10^4)
    xticks=get(gca,'XTick');
    set(gca,'XTickLabel',xticks/1000)
    
    %     colormap(turbo(18))
    %     colormap(cmap)
%     if aa==6
        cbh=colorbar;
        cpos=get(cbh,'pos');
        set(cbh,'pos',[cpos(1)+.04 cpos(2)+.2 cpos(3) .5*cpos(4)])
        ylabel(cbh,'dbZ','fontsize',18,'fontweight','bold')
%     end
    hold on;
    %     ah=area(xvnew,squeeze(terrainint(yvnew==0,:)));
%     ah2=area(xvnew-10^4,nanmax(terrainint(abs(yvnew)<5000,:),[],1));%,'-k','linewidth',2);
% %     ah=area(xvnew,nanmean(terrainint(abs(yvnew)<2500,:),1));%,'-k');
%     ah2.FaceColor=[.8 .8 .8];
%     ah2.LineWidth=3;
%     ah2.FaceAlpha=.5;
%     %         plot(xvnew,nanmin(terrainint(abs(yvnew)<5000,:),[],1),'-k','linewidth',2);

    ah=area(xvnew-(-2.8*10^4),nanmean(terrainint(abs(yvnew)<2500,:),1));%,'-k');
    ah.FaceColor=[.7 .7 .7];
    ah.LineWidth=3;
    daspect([1 .4 1]);
%     xlim([-30 20]*10^3)
    ylim([00 12000])

%     ah3=area(xvnew-(-1.5*10^4),nanmin(terrainint(abs(yvnew)<5000,:),[],1));%,'-k','linewidth',2);
%     ah=area(xvnew,nanmean(terrainint(abs(yvnew)<2500,:),1));%,'-k');
%     ah3.FaceColor=[.8 .8 .8];
%     ah3.LineWidth=3;
%     ah3.FaceAlpha=.5;
    %fire perimeter data
    [val,ttidx]=min(abs(ftimer-timerot(ttt)));%datenum(2020,9,9,1,10,0)));

    %note these data are in x,y coordinates relative to KBBX, whereas the other
    %radar data is relaive to KDAX. Thus we need to add an xoffset and yoffset
    %to these data
    StationLatitudekbbx              = 39.4961;
    StationLongitudekbbx             = -121.6317;
    StationElevationInMeterskbbx    = 53;
    xoff=-(StationLongitude-StationLongitudekbbx).*scalefactor;

    xnowl=(fire_cloud(ttidx).xperim./scalefactor)+StationLongitudekbbx;
    ynowl=(fire_cloud(ttidx).yperim./111180)+StationLatitudekbbx;

    xnow=(xnowl-StationLongitude).*scalefactor;
    ynow=(ynowl-StationLatitude).*111180;
    xrot=((xnow-x0).*cosd(theta)) + ((ynow-y0).*sind(theta));
    yrot=-(xnow-x0).*sind(theta) + ((ynow-y0).*cosd(theta));
    %     fz=F2(xrot,yrot);
    fz2=F2(xrot,zeros(size(xrot)));

    tcross=squeeze(terrainint(yvnew==0,:));
    tinner=inpolygon(xvnew,zeros(size(xvnew)),xrot,yrot);
    sp3=scatter(xrot-(-2.8*10^4),fz2,150,'sr','filled')
    scatter(xvnew(tinner)-(-2.8*10^4),tcross(tinner),150,'sr','filled')
    grid on; box on;
    set(gca,'layer','top','linewidth',2,'fontsize',18,'fontweight','bold')
%     cbh=colorbar;

    xlabel('Dist [m]')
    ylabel('Alt [m]')

    grid on; box on; set(gca,'layer','top','linewidth',2,'fontsize',18,'fontweight','bold')
    hold on
%         text(-1000,11000,strcat(datestr(timerot((ttt)),'HHMM'))),'-',datestr(timer(tidx(ttt+2)),'HHMM'),{' UTC'}),'fontsize',18,'fontweight','bold')
    pause(2);
title(datestr(timerot(ttt)))

    %% now add cross section line to the plan-view map
    %     subplot(1,6,5:6);cla; hold on;
    axx=axes('Position',[.59 .55 .28 .25]); hold on;
  xlim([-121.5 -120.8]);
    ylim([39.5 39.8])
    plot_google_map('MapType','terrain');
    hold on;
    xxl=(squeeze(XX(:,:,1))./scalefactor)+StationLongitude;
    yyl=(squeeze(YY(:,:,1)./111180))+StationLatitude;
    ztidx=find(timer==timerot(ttt))
    dbzmap=squeeze(nanmax(dbzmaster(ztidx,:,:,1:5),[],4));
    dbzmap(dbzmap<10)=nan

    pp=pcolor(xxl,yyl,(dbzmap)); shading flat;caxis([10 45])
    
    pp.FaceAlpha=.75;
    hold on;
    %           contour(squeeze(xxl(:,:,1)),squeeze(yyl(:,:,1)),medfilt2(dbzmap,[5 5]),[10:2:40],'linewidth',2);
    %     contour(xterrain(1:1:end,1:1:end),yterrain(1:1:end,1:1:end),zall(1:1:end,1:1:end),[0:200:3000],'color',[.7 .7 .7]);
    hold on;
    colormap(BEACHIE_SMOKE)

    xcrosser=unique(xxx)
    ycrosser=zeros(size(xcrosser));
    xorg=x0+xcrosser.*cosd(theta)-ycrosser.*sind(theta);
    yorg=y0+ycrosser.*cosd(theta)+xcrosser.*sind(theta);
    hold on;
    scatter(xnowl,ynowl,'*r')

    plot((xorg./scalefactor)+StationLongitude,(yorg./111180)+StationLatitude,'--k','linewidth',3)

    daspect([1 .78 1])
%     cbh=colorbar; ylabel(cbh,'dbZ')
%     title(strcat(datestr(timer(ztidx)-7/24,'HH:MM PM')));%,{' - '},datestr(timer(ff)-7/24,'HH:MM PM'),'{ PDT}'))

    hold on;
    xlim([-121.5 -120.8]);
    ylim([39.5 39.85])
    grid on; box on;set(gca,'layer','top','linewidth',2,'fontsize',15,'fontweight','bold')
    pause(3);
    return
    ofile=strcat('bear_cross_',datestr(timerot(ttt),'yyyymmddHHMM'),'.png')
    export_fig(sprintf('%s',ofile),'-m2.5');

end
return